Controlling the energy flow in nonlinear lattices: a model for a thermal rectifier. 
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We address the problem of heat conduction in 1-D nonlinear chains; we show that, acting on 
the parameter which controls the strength of the on site potential inside a segment of the chain, we 
induce a transition from conducting to insulating behavior in the whole system. Quite remarkably, 
the same transition can be observed by increasing the temperatures of the thermal baths at both 
ends of the chain by the same amount. The control of heat conduction by nonlinearity opens the 
possibility to propose new devices such as a thermal rectifier. 



In recent years a renewed attention has been directed 
to the energy transport in dynamical systems, a prob- 
lem which has been denoted by Peierls as one of the 
outstanding unsolved problems of modern physics [Q. 
These efforts mainly focused on the possibility to derive 
the Fourier law of heat conduction on purely dynamical 
grounds without recourse to any statistical assumption 
|p[-pT[. In spite of relevant progresses, several problems 
remain open and we are still far from a complete under- 
standing nii^i- 

In this paper we investigate a different and impor- 
tant problem namely the possibility to control the en- 
ergy transport inside a nonlinear ID chain connecting 
two thermostats at different temperatures. We show that 
we can parametrically control the heat flux through the 
system by acting on a small central part of the chain. 
Even more interestingly, we show that it is possible to 
adjust the heat flux by varying the temperatures of both 
thermostats, keeping constant the temperature differ- 
ence. Thus we provide a simple mechanism to change 
the properties of the system, from a normal conductor 
obeying Fourier law, down to an almost perfect insula- 
tor. Controlling heat conduction by nonlinearity opens 
new possibilities, such as the design of a lattice that car- 
ries heats preferentially in one direction, i.e. a thermal 
rectifier. 

We consider the Hamiltonian 



H 



E ^ + VM + \K{yn-yn^if (1) 

n=l,N 



which describes a chain of N particles with harmonic 
coupling of constant K and a Morse on-site potential 
Ki(yri) ~ Z3„(e~""*'" — 1)^. This model was introduced 
for DNA chains where m is the reduced mass of a base 
pair, yn denotes the stretching from equilibrium position 
of the hydrogen bonds connecting the two bases of the 
71-th pair and p„ is its momentum Jl2|-|l^ . In the context 
of the present study, model (Q) can simply be viewed as a 
generic system of anharmonic coupled oscillators, the on- 
site potential arising from interactions with other parts 



of the system, not included in the model. The Morse po- 
tential is simply an example of a highly anharmonic soft 
potential, which has a frequency that decreases drasti- 
cally when the amplitude of the motion increases. 

In this paper we consider the out-of-equilibrium prop- 
erties of model (|]) by numerically simulating the dynam- 
ics of the N particle chain, coupled, at the two ends, with 
thermal baths at different temperatures Ti and T2. We 
thermalize at Ti and T2 the first and the last L parti- 
cles by using Nose-Hoover thermostats chains 16 



a Langevin description when we investigate cases very 
far from equilibrium. We then integrate the differential 
equations of motion both for the thermostats and for the 
chain with a 4th order Runge-Kutta method as described 
in 11^. The baths temperatures Ti, T2 are never large 
enough to drive the system beyond the thermal denatu- 
ration temperature Tc above which the mean value of ?/„ 
diverges |Q. 

We compute the temperature profile inside the sys- 
tem, i.e. the local temperature at site n defined as 
T„ — m{yn^), where ( ) stands for temporal average, 
and the local heat flux J„ — K{yn{yn — yn-i)) The 
simulations are performed long enough to allow the sys- 
tem to reach a non-equilibrium steady state with local 
thermal equilibrium where the local heat flux is constant 
along the chain. 

As a preliminary step we have considered the homo- 
geneous case in which _D„ = Z?, a„ = a, n = 1, ..N . 
Here, as expected, we have detected a temperature gra- 
dient inside the chain, as shown in Fig. |l| (circles), and 
we have verified that the system obeys the Fourier law of 
heat conduction. For fixed temperatures Ti and T2, and 
a chain length varying from = 64 to iV = 1024 the 
heat flux evolves as 1/A^ as expected for a system with a 
well defined, finite, thermal conductivity. 

We now divide the chain between the thermostats in 
three regions in which D„ takes different values. In 
the left and right regions, n = 1, . . . (A — M)/2 and 
n = (A^ -I- M)/2 -I- 1, . . . , A^, respectively, we set £)„ = £), 
while in the M sites of the central region D„ = Di; 
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an = OL for the whole chain. Typically, in our simula- 
tions, we take D = 0.5, a = 1, K = 0.3, m = 1, N = 128, 
L = 16 and M = 8 while Di has been varied from 0.5 
to 1.2. As shown by the numerically computed temper- 
ature profiles of Fig. the results are now strikingly 
different. The temperature profile changes as Di — D is 
increased until, for large enough Di, the central region 
behaves as an insulator, the left and right regions being 
thermalized at Ti and T2 respectively, and the heat flux 
drops to negligible values ~ 10~^. Therefore the change 
in Di can induce a conductor-insulator transition con- 
firmed by Fig. ^ which shows that the averaged heat flux 
J decreases by two orders of magnitude when Di — D 
increases from to 0.7. 
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FIG. 1. Temperature profile of the chain with parameters 
D = 0.5, N = 128, 1, K ^ 0.3, M = 8, and Ti = 0.16, 
T2 — 0.15. The Morse potential constants are Di — D — 0.5 
- homogeneous system- (circles), Di = 1.2 (diamonds) and 
Di = 0.8 {+) 

In order to understand this phenomenon, let us con- 
sider a linearized version of model (1), which is obtained 
by linearizing the Morse potential around j/„ = 0. The 
Hamiltonian now writes: 

^= E ^+Dnyl + \Kiyn-yn-i? (2) 

i=7l,N 

where = DnCx^ and _D„, a„ are defined above. 

Numerical simulations with Hamiltonian (||), show 
that this system displays the same conductor-insulator 
behavior in a clearer and sharper way. Indeed, due to 
integrability, there are some typical pathologies, such as 
the presence of plateaus instead of gradients in the tem- 
perature profile, as predicted by |l^. This is why the 
transition from conducting to insulating regime is even 
sharper than in the nonlinear case (Fig. |^ (circles) ) until 
a cut-off is reached. 

The equations of motion for the system (|^), without 
thermostats, have plane wave solutions, ?/„(i) = e*'^"~*"*, 
where the frequency uj and momentum k satisfy the dis- 
persion relation lo^ — 2K + 2D — 2K cos k so that the 
allowed frequencies are included in the so called phonon 
band 2D <uj^ <2D + AK. 
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FIG. 2. Heat flux as a function of parameter Di — D; the 
squares refer to the Morse potential, while the circles refer to 
the corresponding harmonic limit. Here A'^ = 128, D = 0.5, 
K = 0.3, a = 1, L = 16 

Let us now consider the harmonic and inhomogeneous 
chain, with Di ^ D in the central region. In order 
to propagate through the system a wave should satisfy 
both dispersion relations uj"^ = 2K + 2D — 2K cos k and 
= 2if + 2Di - 2/i:cosA:i. U Di D only a fraction 
of the allowed frequencies of the waves coming from the 
left or right regions can satisfy both dispersion equations 
with real fci, and therefore also propagate through the 
central region. For the remaining frequencies, the second 
equation can only be solved by imaginary fci = iki which 
gives rise to states ?/„ ~ g-fcm which decrease exponen- 
tially from the edges of the central region and therefore 
cannot propagate through it so that their contribution to 
the total heat flux J is exponentially small. By increas- 
ing the separation of the two phonon bands Di — D, the 
fraction of frequencies which can give rise to plane waves 
inside the entire system falls down and reaches zero as 
Di — D > 2K, where the band separation is larger than 
the band width, and the bands do not overlap anymore. 
In this situation, waves coming from the thermostats are 
reflected back from the central region, the heat flux is ex- 
ponentially small and the system acts as an almost per- 
fect insulator. In the case of partially overlapping bands 
one obviously expects an intermediate situation. The 
dependence of the heat flux on the band overlapping is 
shown by the lower curve in Fig. |^ and its variation can 
be calculated with the above analysis. 

The addition of a nonlinear term to the Hamiltonian 
does not change the substance of the above argument. 
The numerically observed smoother decay of the heat 
flux vs. Di — D given by the upper curve in Fig. ^ is due 
to the fact that plane waves are not exact solutions for 
the nonlinear case and the presence of solitary excitations 
like breathers may play an important role. 

It is now interesting to investigate whether one can 
control the transport properties of the chain by acting on 
the temperature of the heat baths without modifying the 
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parameters of the interaction or of the on site potential. 
To this end we recaU that, at constant T, the nonhnear 
model can be approximated by a fully harmonic Hamil- 
tonian, with temperature dependent parameters fl^{T), 

<f(r) m 



not obtained and therefore the system does not become 
a perfect conductor. 



n=l,N 



Vn 



(3) 



An effective-phonon analysis, shows that, while 
$(T) = if is temperature independent, ^{T) decreases 
as T increases The lowering of the band vs. the 

temperature is of course determined by the parameters 
of the system, a„, D„. Typical band profiles are shown 
in the inset of Fig. |^, for different parameter choices. 




FIG. 3. Temperature profiles for system (|l|) at three differ- 
ent values of T2 for the same temperature gradient T2 — T\. 
Parameters are iV = 128, M = 8, ai = 2, a = 1, Di = 0.375, 
D = 0.5; Ta = 0.09 ( circles), Ta = 0.19 (x) and T2 = 0.29 
(diamonds). 

Inset: Phonon band profiles vs. temperature; the continuous 
lines represent the phonon band for a — 2, D = 0.375; the 
dashed lines D = 0.5, a — 1. The full phonon band is shown 
with bandwidth 4K = 1.2 

We may now choose parameters a = 1, D = 0.5 
(dashed lines in the inset of Fig. (||)) for the left and 
right regions, and ai = 2, Di = 0.375 (fidl line in the 
inset of Fig. (^) for the central region. From the inset of 
Fig. (||) it is clearly seen that, at low T, the phonon band 
of the central region does not overlap with the band of 
the other two side regions. In this situation, the heat flux 
is nearly zero and the temperature profile shows a clear 
insulating behavior, with first region thermalized at tem- 
perature Ti and third region at temperature T2 as shown 
in Fig. 1^ (circles). As we increase both temperatures Ti 
and T2, keeping \Ti — constant, the bands start to 
overlap, the temperature profile exhibits a gradient and 
the heat flux increases drastically (Fig. |^). 

As temperature increases, nonlinearity becomes more 
and more important and the effective phonon analysis be- 
comes less and less meaningful. Nevertheless, numerical 
analysis shows that the above argument of band over- 
lap still remains valid even if a clean thermal gradient is 
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FIG. 4. Schematic picture of the phonon bands in a "ther- 
mal rectifier" for two directions of the temperature gradient. 
The bands in the left and right weakly anharmonic regions 
do not change significantly with the orientation of the gradi- 
ent. In the central part, when the high temperature side is on 
the right, the band evolves in space as shown by the shaded 
region, while, if the high temperature side is on the left, the 
band evolves according to the dashed lines. 

Once the general mechanism of the thermal conduction 
in a composite nonlinear lattice has been understood, this 
opens many possibilities. For instance one can design a 
thermal rectifier, as schematized on Fig. ^. A strongly 
nonlinear region is sandwiched between two weakly an- 
harmonic left and right domains. In the presence of a 
thermal gradient in the central part, the effective phonon 
frequencies evolve in space in a way that depends on the 
orientation of the gradient. This can provide either a 
good matching of the bands at the interfaces, leading 
to a thermal conduction across the system, or a com- 
plete mismatch leading to poor conduction. This picture 
is oversimplified for various reasons: the self-consistent 
phonon method is approximate, there is not a homoge- 
neous gradient in the system due imperfect thermal con- 
tact at the interfaces, and, in spite of nonlinear mixing, 
in the conducting direction, the frequencies injected at 
one side of the central region do not shift gradually, fol- 
lowing the phonon band, to match the band of the other 
outer region. 

In spite of these restrictions Fig. H shows that it is pos- 
sible to build a simple thermal rectifier where the thermal 
fiux changes by a factor of about 2 when the direction of 
the thermal gradient is reversed. In fact thermal rectifiers 
are already known: a layer of fluid, heated from below 
or from above yields very different thermal flows because 
the former case leads to convection and the latter does 
not. To our knowledge, it is nevertheless the first time 
that a "solid system" (a highly simplified model of a solid 
indeed) shows such a behavior. 

In this work we have studied a specific nonlinear lat- 
tice but the results rely on ideas of general validity, the 
matching of phonon bands or the frequency shift associ- 
ated to nonlinearity, so that one can expect that various 
physical systems could exhibit similar behaviors. The 
simple model that we have chosen to introduce these 
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ideas involves an on-site potential which describes the 
effect of an external substrate not explicitely included in 
our calculations. In an actual device, such a substrate 
would introduce a second channel for thermal energy 
transport, operating in parallel with the channel stud- 
ied here. This could reduce the efficiency of the rectifier, 
unless the thermal conductivity of the substrate is signif- 
icantly lower than the conductivity of the lattice of non- 
linear oscillators. Any solid thermal insulator or weak 
conductor, can, in principle, be used as a substrate. A 
prototype could be a disordered harmonic lattice pl|]. 



open problem. 
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FIG. 5. Temperature profiles in a "thermal rectifier" for 
two opposite orientations of the thermal gradient. The 
dash-dotted lines show the borders of the different regions 
in the lattice. The thermostats have a size L = 8, the central 
region a size AI — 128, and the left and right regions con- 
tain 64 particles. The coupling constant is K — 0.18, and 
the parameters of the Morse potential are D = 4.5, a — 0.5, 
Di = 0.7, ai = 1.4, D' = 2.8, a' = 0.5 in the left, cen- 
tral and right regions respectively. The temperatures of the 
thermostats are T = 0.1 and 0.45. When the high temper- 
ature is on the right side of the lattice, the average flux is 
J — 0.146 10"'^ while, when the thermal gradient is reversed 
the flux drops to J = 0.755 10~*. Note in this case the dis- 
continuities at the interfaces between region which attest of 
the bad energy transfer at these points. 

Among the possible applications of this study, one can 
think of biological molecules. Controlling the energy flow 
in biomolecules is important, for instance when the en- 
ergy of ATP hydrolysis is released locally in a molecular 
motor to be used elsewhere after some delay ||2^]. Elec- 
tronic states could be involved, but we show here that a 
control of the flow of thermal energy can also be achieved 
through the strong nonlinearities which exist in biolog- 
ical molecules and could even be tuned by events such 
as a conformational change. Experiments show that the 
energy flow in proteins can be extremely slow pl[|, pro- 
viding a weakly conducting substrate. On the other hand 
our results show that a controlable conducting chan- 
nel attached to this substrate is conceivable. Whether 
biomolecules exploit this possibility is still a completely 
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